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1. Introduction 

Understanding time evolution of an open quantum system of many interacting particles 
is of primary importance in fundamental problems of quantum physics, such as 
decoherence [HE] and closely related quantum measurement problem 0,1!], quantum 
computation [5j [6], or the problem of computation of non- equilibrium steady states 
(NESS) in quantum statistical mechanics [3 EJ [9]. Even though application of the 
methods of Hamiltonian dynamical systems and ergodic theory to quantum systems 
out of equilibrium gives many promising results JTUJ HH [12] , the field of open quantum 
systems is still lacking non-trivial explicitly solvable models, as compared to studies 
of closed (isolated) quantum systems where we know a large body of the so-called 
completely integrable systems [131 [H]. Examples of explicitly solvable models of master 
equations for open quantum systems are limited to quite restricted models of a single 
particle, single spin or harmonic oscillators (see e.g. [T5 | IT6l \TT\). 

In this paper we show that the generator of the master equation of a general 
quadratic system of n interacting fermions which are coupled to a general set of 
Markovian baths, specified in terms of Lindblad operators which are linear in the 
fermionic variables - the so called quantum Liouville super-operator (or Liouvillean) - 
can be explicitly diagonalized in terms of 2n normal master modes, i.e. anticommuting 
super-operators which act on the Fock space of density operators. This can be 
understood as a complex (non-canonical) version of the Bogoliubov transformation [18] 
lifted on the operator space, and has very powerful consequences: (i) The NESS of the 
master equation can be understood as the 'ground state' normal mode of the Liouvillean, 
whereas the long time relaxation rate is given by the eigenmode closest to the real 
axis, (ii) The covariance matrix of NESS can be computed explicitly in terms of the 
eigenvectors of 4n x 4n antisymmetric complex matrix. It can be used to completely 
express physical observables in NESS, such as particle/spin densities, currents, etc. 

We demonstrate the power of this novel method by applying it to the problem of 
heat and spin transport far from equilibrium in nearest neighbor Heisenberg XY spin 1/2 
chains subject to a transverse magnetic field. As a result we reproduce ballistic transport 
in the integrable spatially homogeneous case (see e.g. [HI [20l [211 [221 [231 [2H 12] for 
related recent studies of quantum thermal conductivity in one dimension), and predict 
ideally insulating behaviour (at all temperatures) in a disordered case of spatially random 
interactions/field. Apart from obtaining numerical results which go by far beyond 
what was so far accessible by direct numerical solution of the many-particle Lindblad 
equation, either directly or by means of quantum trajectories |15j . we also obtain two 
notable analytical results in the spatially homogeneous (non-disordered) case: (i) We 
compute the spectral gap of the Liouvillean i.e. the rate of of relaxation to the NESS 
and show that it scales with the inverse cube of the chain length, (ii) We construct 
evanescent normal master modes of the Liouvillean, for long chains, by which we explain 
quantitatively the exponential falloff of energy density or temperature profiles near the 
bath contacts. 
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The paper is organized as follows. In section [2] we shall outline a general method 
for the diagonalization of the Liouvillean super-operator for finite quadratic open Fermi 
systems and an explicit construction of NESS. In section [3] we illustrate the method 
by working out a simple example of a single fermion or a two level quantum system in 
a bath. In section H] we demonstrate the usefulness of the new method by applying it 
to quantum transport in XY spin chains. In section [5] we discuss possible alternative 
applications and generalizations of the method and reach some conclusions. 



2. General method of solution for the Lindblad equation 

The general master equation governing time evolution of the density matrix p(t) of 
an open quantum system, preserving trace and positivity of p, can be written in the 
Lindblad form [25l [17] as (we set h = 1) 

f t = tp := -i[H, p] + J2 ( 2L »P L l ~ HW, P}) (1) 

n 

where H is a Hermitian operator (Hamiltonian), [x,y] := xy — yx, {x,y} := xy + yx, 
and are arbitrary operators representing couplings to different baths (at possibly 
different values of thermodynamic potentials). We are now going to describe a general 
method of explicit solution of ([I]) for a quadratic system of n fermions (or spins 1/2) 
with linear bath operators 

2n 

H = WjHj k Wk = w-Hw (2) 
j,k=i 

2n 

L »=Y1 l ^ w i ~ -m " — ( 3 ) 

3=1 

where Wj, j = 1,2, .. . , 2n, are abstract Hermitian Majorana operators satisfying the 
anti-commutation relations 

{vjj,w k } = 26j lk j,k = l,2,...,2n (4) 

and generate a Clifford algebra. Thus, 2n x 2n matrix H can be chosen to be 
antisymmetric H T = — H. Throughout this paper x = (x\, X2, ■ ■ -) T will designate a 
vector (column) of appropriate scalar valued or operator valued symbols x k - 

Two notable examples, to which our formalism is immediately applicable, are: (i) 
canonical fermions c m , m = 1, 2, . . . , n, 

w 2m -i = c m + c4 w 2m = i(c m - 4J (5) 
or (ii) spins 1/2 with canonical Pauli operators a m , m — 1, 2, . . . , n, 

w 2m -i = < n = a ™ n ( 6 ) 

m' <m m' <m 

Here we are not concerned with physical criteria for the validity of the so-called 
Markovian approximation under which eq. (pQ) is derived, so we shall make no 
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assumptions on the smallness of the bath coupling constants l^j. We merely consider the 
Lindblad equation ([!]) as a possible parametrization of an important subset of Markovian 
completely positive quantum channels and demonstrate its complete solvability for 
quadratic systems. Note that generalization of our formalism to explicitly time dependent 
Hamiltonians H(t) and Lindblad operators L^t), generating more general and possibly 
non-Markovian open system dynamics, is straightforward. See e.g. [26] for a discussion 
of Markovianity. 



2.1. Fock space of operators 

We begin by associating a Hilbert space structure x — > \x) to a linear 2 2n = 4™ 
dimensional space /C of operators, with a canonical basis |P„) with 

P ai ,a 2 ,...,a 2n := w?w? ■ -.w%r a, e {0, 1} (7) 

orthonormal with respect to an inner product 

(x\y) = 2~ n Wy (8) 

The form of the canonical basis of the operator space K, suggests that it is just a usual 
Fock space with an unusual physical interpretation. Namely we can define the following 
set of In adjoint annihilation linear maps dj over /C 

CjlPa) = S^^WjPa) (9) 

and derive the actions of their Hermitian adjoints - the creation linear maps &, 

(Ptf\c]\Pa) = (Pa I C? I -Pa')* = <*a$,l ( P » \ w j P a')* = 8 aj ,o( P a'\WjP^): 

c]\Pa) = 5 ajfi \w 3 Pa) (10) 

Straightforward inspection then shows that they satisfy the canonical anticommutation 
relations 

{c j ,c k } = {cj,c\} = 5 jik j, k = 1,2, . . .,2n (11) 

The key is now to realize that the quantum Liouville map £ defined by eqs. ( 1112131) can 
be written as a quadratic form in adjoint Fermi maps dj, ct (or for short, a-fermions) . [|| 

First, we consider the unitary part of Liouvillean 

C Q p:=-i[H,p] (12) 

Since K, is a Lie algebra, one defines the adjoint representation of a Lie derivative for 
an arbitrary A G K, back on JC as, &dA\B) := \[A, B]). It is now straightforward to 

| Throughout this paper Dirac's bra-ket notation shall be used only for a Hilbert space K. of physical 
operators, including density operators, in a sense of GNS construction although here all spaces will 
be finite dimensional. Symbols with a hat shall designate linear maps over the operator space K,. For 
instance, we note a key distinction between physical fermions c m ([5]) and a-fermions Cj (J5J)- 
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compute the action of a Lie derivative of a product of two Majorana operators on an 
arbitrary basis element 

ad WfWfc | Pa) = \WjWkPa) - \P<y_WjW k ) = 

= 2 (<5 Q . ,1^,0 + S a . }0 5 aktl )\wjW k Pa) = 

= 2(C,4 + c]c,)|Pa) = 2(c]6 k - 4c,-)|Pa) (13) 

Extending this relation by linearity to an arbitrary element of /C, it follows that for an 
arbitrary quadratic Hamiltonian ([2]) its Lie derivative has a very similar quadratic form 
in a-Fermi maps 

2n 

Co = -iad# = -4i c)H jk c k = -4ic f • He (14) 

j,k=i 

It is worth stressing here that for an arbitrary (complex) matrix H, £ (fl4l) conserves 
the total number of a-fermions J\f := ^ ■ c^Cj — $ ■ c, namely [Co, J\f] = 0. 

Second, we consider the action of the Lindblad maps 

2n 

C^p := 1L^pL\ - {L\ L L^ p} =Y1 l M l l,A,kP (15) 

j,k=i 

where we write Cj tk p '■= 2wjpw k — w k Wjp — pw k Wj. Again we proceed by computing 
the actions of Cj tk on elements of the canonical basis of operator Fock space /C. In 
order to do so, it is crucial to observe that the question whether Wj commutes or 
anticommutes with depends on the number of a-fermions \a\ := YmcLi a k i n \Pa)i 
namely PgWj = (— l)'-' +Qj WjPa, and hence 

Cj, k \Pa) = [2(-lp +a " Wj w k - w k w 3 - {-l)^ +ak w kWj \ | P„) (16) 



Observing that 



\wjPa) = (cJ + c^lPa) (IT) 

{-l) a >\w j Pj = $-c j )\Pj (18) 

(-l)N|Pa) = exp(i7rA0lP^ (19) 
one derives from ({TBI) the general expression for Cj tk 

A'.fe = (l + exp(i7rA0 ) (%c]c\ - d]c k - c\c^j 

+ (l - exp(iTrAO) (2cjC k - cfc\ - c k cf\ (20) 

Obviously, the maps Cj tk , and hence also the total Lindblad part of Liouvillean Y^u^-vi 
do not conserve the number of a-fermions. But they conserve its parity i.e. the product 
of any two creation/annihilation a-Fermi maps commutes with the parity operation 
V = exp(mAf), with respect to which the operator space can be decomposed into a 
direct sum /C = /C + © /C~, and even/odd operator spaces are orthogonally projected as 
IC^ = (i ±exp(i7rA/"))/C. Thus the positive parity subspace /C + is a linear space spanned 
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by | -P a ) with even \a\. All the maps Cj tk now act separately on /C 1 * 1 , Cj^K.^ C /C ± . For 
example, the maps defined on even parity subspace are indeed quadratic in a-fermions 

4fek+ = 4c}4 - 2cjc fc - 2c\cj (21) 

In this paper we shall focus on physical observables which are products of an even 
number of Majorana fermions - operators Wj - so we shall in the following discuss 
only Liouville dynamics on the subspace /C + . The extension to the dynamics of odd 
observables should be straightforward. 

Putting the results ( TT2lll4l)l 5|f2Tj) together we arrive at the final compact quadratic 
representation of the Liouvillean C + := £|ye+ 

t + = -2 c f ■ (2iH + M + M T ) c + 2 c f ■ (M - M T ) c f (22) 
where M is a complex Hermitian matrix parametrizing the Lindblad operators 

M jk = Y, l Ak (23) 

2.2. Reduction to normal master modes 

Next we want to show that the representation ( 1221) allows us to reduce it further by 
a linear transformation of the set of maps {cj,c[;j = 1,2, ...,2n} to normal master 
modes (NMM) in terms of which the complete spectrum of the Liouvillean, as well as 
its eigenvectors, can be explicitly constructed; in particular the zero-mode eigenvector 
which is just the physically relevant NESS. 

In fact we proceed in analogy to Lieb et al. [H] and define An adjoint Hermitian 
Majorana maps d r = aj, r = 1, 2, . . . , 4n: 

1 i 

a 2j _! = -j=(cj + ct) a 2j = -j=(cj - ct) (24) 

satisfying the anti-commutation relations 

{a r , a s } = 5 r>s (25) 
in terms of which the Liouvillean (1221) can be rewritten as 

C + = a-Aa-A Q l (26) 
where A is an antisymmetric complex An x An matrix with entries 



A-2j-l,2k-l — 


- 2iH jk - M jk + M kj 




A 2 j-l,2k = 


2iM kj 




A2j,2k-1 = 


- 2iM jk 




A-2j,2k = 


- 2iH jk + M jk - M kj 


(27) 



1 is an identity map over /C and Aq is a scalar 



2/i 



A = 2 Mjj = 2 tr M (28) 

3=1 
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Obviously, the bi-linear Liouvillean (1261) cannot be brought to a normal form with a 
linear canonical transformation since the matrix A - which shall in the following be 
referred to as a shape matrix of Liouvillean - is not anti-Hermitian like in Hamiltonian 
systems. So we should proceed in more general terms. 

We first recall few facts about complex antisymmetric matrices of even dimension. 
If v is a right eigenvector Av = (3y_ with complex eigenvalue (3, then y_ is also a left 
eigenvector with eigenvalue — j3, A T y_ = —Av_ = —j3v. Hence eigenvalues always 
come in pairs (3, —f3. Let as assume that A can be diagonalized%, i.e. there exist 
4n linearly independent vectors v r ,r = l,...,4n with the corresponding eigenvalues 

Pli ~Pl, P2, —02, • • • , Pint — Pin, 

Avjy-i = tei-i A ^2j = -PjV 2j (29) 
ordered such that Re Pi > Re P2 > ■ ■ ■ > Re /?2n > 0. The 2n complex numbers Pj shall 
be referred to as rapidities. It is easy to check that we can always choose and normalize 



v r such tha 

1 

v r 'V-s = Jrs where J := a x <S> l^n — 



V: 

Let V be 4n x 4n matrix whose rth row is given by v_ r , V rs := v, 
rewrite as 

tT 



1 











1 






1 





(30) 



\ 



/ 



Then eqs. ( 1291130 



AV J 
VV T 



V T D 



where D :- 



diag{/3i, -fa, /3 2 , -fa, • • • , fizn, -fan} (31) 

(32) 

Expressing V T in terms of (132]) and plugging the result into eq. (13T!) we arrive at a very 
convenient canonical form of a generic complex antisymmetric matrix A 



A = V T AV where A = DJ 



( 


Pi 








-Pi 




























-Pi 






V 



\ 



(33) 



Now we apply decomposition ( |33l) to the Liouvillean ( |26i) 
C + = a ■ V T AVa - A 1 = (Va) • A(Va) - A 1 
Let us define the NMM maps b := (bi,b[, b 2 , b' 2 , . . . , 6 2 n, b' 2n ) : = 



Va or 



v 2j -i ■ a 



v.2j ■ a 



(34) 



(35) 



§ It is not known at present whether explict form P?)) guarantees diagonalizability of any such A. Note 
that one can construct certain types of complex antisymmetric matrices with degenerate eigenvalues 
which cannot be diagonalized [27] , 

|| For a non-degenerate rapidity spectrum {f3j} the proof of this statement is a trivial consequence of 
antisymmetry A = — A T , whereas in case of degeneracies it can be shown that one can always choose 
appropriate linear combinations of eigenvectors. 
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We note that due to (I25|I30I) NMM maps satisfy almost- canonical anti-commutation 
relations 

{b J ,b k } = {hX} = 5 hk {b' 1 ,b' k } = (36) 

namely bj could be interpreted as annihilation map and as a creation map of jth 
NMM, but we should note that b'j is in general not the Hermitian adjoint of bj [28J. In 
terms of NMM the Liouvillean (l34"j) now achieves a very convenient normal form 

2n 

£+ = -2^/3^4-501 (37) 

where Bq = Aq — Yl^iPj- We shall later show that the constant Bq is in fact equal to 
0. 



2.3. N on- equilibrium steady states and a complete spectrum of the Liouvillean 

The Liouvillean can always be represented in terms of a large but finite 4 n x 4 n matrix. 
We shall now outline the procedure of complete construction of its spectrum in terms 
of NMM which are easy to calculate in terms of diagonalization of 4n x 4n matrix A as 
described in the previous subsection. 

We proceed by constructing the Liouvillean 'vacuum'. From the representation 
(J22l) it follows immediately that (1| = (P(o,o...,o)| is left-annihilated by C + , = 0, 

or equivalently = 0. So we have just shown that is always an eigenvalue of C + , 

hence there should also exist the corresponding right eigenvector |NESS), normalized as 
(1|NESS) = trpNESs = 1) which represents physical NESS, i.e. stationary solutions of 
the Lindblad equation ([I]) 

£+|NESS) = (38) 

Let us define NMM number maps as Afj := b'jbj. From eqs. (1361) it follows that Afj satisfy 
a projection property Aff = A/}, so they are diagonalizable since no nontrivial Jordan 
block could satisfy the projection property. Furthermore, Afj are mutually commuting 
[J\fj,Afk\ = 0, so they can be simultaneously diagonalized and there should exist a 
vacuum state on which all Afj have value 0. It follows from the stability of completely 
positive evolution ([TJ) that all eigenvalues A of t + should obey Re A < 0, and since by 
assumption Re/?j > 0, (1| and |NESS) should be the left and right vacua which are 
simultaneously annihilated by NMM maps 

= ^|NESS) = (39) 

and hence also A/}|NESS) = 0. Thus we have also shown that the NMM representation 
(1371) is only consistent if -Bo = so we find an interesting sum rule for rapidities 

2n 

J2f3j = 2 tiM (40) 
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The complete excitation spectrum and the corresponding left/right eigenvectors of 
the Liouvillean are given in terms of a sequence of 2n binary integers (NMM occupation 
numbers) v = (y u v 2 , . . . , u 2n ), Vj e {0, 1}, 

(& L U \C + = X E (Q L J £ + |6*) = A„|9*) (41) 



namely 



2n 



K = ~2J2^ ( 42 ) 



<ey = <i|C ■■■K 2 |ej> = &r bp ■ ■ ■ jciness) (43) 

where by construction, left and right eigenvectors satisfy the bi-orthonormality relation 

(e^| e*} = 

TTie mam general results: uniqueness of NESS, rate of relaxation to NESS, and 
expectation values of physical observables 

Given a physical observable X G /C + and an arbitrary initial state with a density operator 
Po £ the time dependent expectation value of X can be written in terms of the 
spectral resolution of the Liouvillean, 

exp(tC + ) =J2^MtK)\&l)(& L J (44) 

namely 

(X(t)) = trXp(t) = tr[xexp(iC + ) Po ] = ^exp(tA i ,)(e^|p X|e2) (45) 

v_ 

We remind the reader that £ + correctly represents physical Liouvillean only on the 
subspace /C + of operators with even number of a-fermions. However, since the dynamics 
is closed on /C + and test physical observable X also belongs to /C + it follows that the 
component of p from fC~ does not contribute to the expectation value (X(t)}. 

Given the exact and explicit constructions developed in this section we can now 
make the following rigorous and useful conclusions, assuming throughout that Liouvil- 
lean shape matrix A is diagonalizable: 

Theorem 1: NESS of Lindblad equation is unique if and only if the rapidity spec- 
trum {/3j} does not contain 0, in our ordering convention, if /3 2n ^ 0. In the opposite 
case, if we have d > 1 vanishing rapidities, then we have a 2 d dimensional convex set of 
non-equilibrium steady states which can be spanned with 0ui 

Theorem 2: An arbitrary initial state with a density operator p G K, converges with 
time to NESS if and only if all rapidities have strictly positive real parts, Re/3j > 0. 
Then, the rate of exponential relaxation to NESS is given by the spectral gap A of the 
Liouvillean which equals A = 2Re/32n- 
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Theorem 3: Assume that the rapidity spectrum does not contain 0, i.e. (3 2n 7^ 0. 
Then the expectation value of any quadratic observable WjWk in a (unique) NESS can 
be explicitly computed as 

(wjW k ) NESS = 8 j>k + (l|c,-c fc |NESS) = (46) 



S 



'j,k 



i 2n 



v 2m,2j-l v 2m-l,2k-l — V2 m ,2j v 2m-l,2k 



m=l 



— 1 V2m,2jV2m-l,2k-l ~ Y ^2m,2j-\^2m-\,2k 



(47) 



The statements of theorems 1 and 2 simply follow from exact and explicit spectral 
decomposition ( T4l2j43ll44l) . 

The proof of theorem 3 is also straightforward: The first expression (|46j) follows 
from the definition of the annihilation maps ([9]) and the explicit representation of the 
density operator of NESS, Pness; m the canonical basis P a . The second, very useful 
equality ( 147|) is then obtained by expressing dj thru f |24l) in terms of NMM maps ( |35l) 
and using the annihilation relations ( 1391) . 

The quadratic correlator of theorem 3 covers many physically interesting 
observables such as densities or currents. However if one needs expectation values 
of more general observables, e.g. an expectation value of a high order monomial 



NESS 



'llc" 1 ^ 2 • • • &2n n |NESS), then one may use a Wick theorem and rewrite 



it as a sum of products of pair-wise contractions fj46l) . 
3. Trivial example: A single fermion in a bath 

In order to illustrate the method and demonstrate convenience of the results derived in 
the previous section we first work out a simple example of a single fermion n = 1 (or 
equivalently, an arbitrary qubit, a two- level quantum system), in a thermal bath. We 
take the most general single fermion Hamiltonian H = —ihwiW2 + const = 2hc^c+ const' 
and the following bath operators (see e.g. [HU 129] ) 



Li = - 

1 2 



Vx(w 1 - iw 2 ) 



L 2 = ^V^iwi + iw 2 ) 



r 2 c t (48) 



where the ratio of coupling constants determine the bath temperature T, r 2 /r! = 
exp(— 2h/T). Leaving out the details of a straightforward calculation, simply following 
the steps of the previous section, we arrive at the following shape matrix of the 
Liouvillean (1261) 

A = -/iR + B r , r _ A) = r + 



where 



R 



/ 1 0\ 

1 

-10 

V o -l o o/ 



Br + ,r_ 



( 




2 L - 


ir 

2 L - 


\ 


-2- r + 





ir 

2 ~ 


ir 

2 - 






-ir 

2 L ~ 





2 1 + 




V-|r_ 


2 L ~ 


-*r + 





/ 



(49) 



(50) 
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and r± := r 2 ± IY Further, we compute NMM rapidities /3i j2 = hT + ± ih and the 
eigenvector matrix 



V 



V 



• r_ 



- + 1 ip i ijr- + i 



(51) 



/ 



Then, using theorem 3 we compute the expectation value of occupation number 



(etc) 



h(wiW 2 ) = r 2 /(r! + r 2 ) which is what we expect in canonical equilibrium. 



4. Non-trivial example: transport in quantum spin chains 



Here we work out a physically more interesting example, namely we construct NESS 
for the magnetic and heat transport of a Heisenberg XY spin 1/2 chain, with arbitrary 
- either homogeneous or positionally dependent (e.g. disordered) - nearest neighbour 
interaction 

n— 1 n 

H =J2 ( J »m+1 + JK^m+l) + Yl h ™ a ™ ( 52 ) 
m=l m=l 

which is coupled to two thermal/magnetic baths at the ends of the chain, generated by 
two pairs of canonical Lindblad operators |29j (similar to (|48j) ) 

l f-r _ . i 



Li 
L 2 



(53) 



2 v z i 2 

where = cr^ ± icr^ and are positive coupling constants related to bath 

temperatures/magnetizations, for example if spins were non-interacting the bath 
temperatures T Lj r would be given with r^/r^' 11 = exp(— 2/i l n /T Li R). 



Applying the inverse of Jordan- Wigner transformation (EJ), cr^ = (— i) m 1 Yl^=i 1 w ji 



-i 



im-l 



(nSi 2 w j )w 2m , 



we rewrite (I52II53P in terms of Major ana fermions 



n-l 



H = - i ^2 (J£ l w 2m w 2m+ i - Jl l W 2 m-lW 2 



m+ 2 ) 



Li = - 
2 



m=l 
L 



i h m w 2m ^ 1 w 2r 

m=l 



(54) 



r^(wi - iw 2 ) L 3 



L 2 



r 2 (wi + iu> 2 ) L 4 



rf (w 2 „_i 



n-l 



itU2n)W 



(55) 



where W = W\W 2 - ■ ■ w 2n is a Casimir operator which commutes with all the elements of 
the Clifford algebra generated by Wj and squares to unity W 2 = 1. Therefore, W does 
not affect the action of bath operators (TT5]) where enter quadratically, so we find 

ti + t 2 = 



rt(^ci + a^c 2 ) 



2iI^c}cJ 



£ 3 + £ 4 



r+(c 2 n-lC2n-l + C 2 „C 2 „) — 2iF^C 



at 



2n-l u 2re 



(56) 
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leading to the bath shape matrix (1501) identical to the single fermion case f)48p . Again, 
carefully following the steps of section [21 we derive the Liouvillean in the form (T261) 
with 4n x An shape matrix, which we write in a block tridiagonal form in terms of 4 x 4 
matrices as 



/Br 



A 



V 



Rf 






Ri 
-h 2 R 

— RT 





R 2 
-h 3 R 








R 



n— 1 







(57) 



and A = r+ + r+, where B L 



rL.R. 
± - = 



r L,R ± pL.R and 



Rr. 



B 



B 



-R^_l Br — /i n R J 

= B r R r a fin terms of fl50l)'). with 



/ 








P 


\ 













P 




- 7 X 











V 





— 7 X 





o / 



(58) 



We are not able to perform a complete diagonalization of the antisymmetric matrix 
A floTl) of the general XY model analytically. For example, even in the spatially 
homogeneous case J£f = J x ' y , h m = h it is not possible to proceed like in the 
classical harmonic oscillator chain where the corresponding matrix is a sum of a 
Toeplitz and a bordered matrix [30J. Namely, in our case A is a sum of a block 
Toeplitz and block bordered matrix and its explicit exact diagonalization remains an 
open problem. However, we stress that even relying on numerical diagonalization of A 
yielding a set of rapidities /3j and properly normalized eigenvector matrix V, represents 
a dramatic progress with respect to previously existing numerical methods which needed 
diagonalization of matrices which were exponentially large in n. We shall later derive 
some exact theoretical and analytical results, explaining results of exact numerical 
computations, in the special case of a homogeneous transverse Ising chain (subsection 
14.11) . and the case of a disordered XY chain (subsection 14.21) for which we shall relate 
NMM to the problem of Anderson localization in one dimension, 

Let us continue by discussing transport observables in the spin chain whose 
expectation values in NESS are easy to calculate. Note that the bulk Hamiltonian 
(|52|) can be written in terms of the two-body energy density operator 



-i W 2m W 2m +l + iJl l W 2 m-lW 2m +2 



ih r 



ih 



-W 2m -iW 2m 



m+l 



W 2m +lW 2m+2 (59) 



H m = 

as H = J2m,Hm- One can derive the local energy current Q m = i[H m ,H m+ i] from the 
continuity equation 

(d/dt)(H m ) = trH m dp/dt = (i[H,H m \) = -(Q m ) + (Q m -i) 

where Q m := i[H m , H m+1 ] 

Qm — 2K'2J m Jm-i-l W 2m-l w 2m+3 + 2J^, Jj, 



(60) 



^+l'U , 2m W 2m+ 4 : — 
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— h m+ iJ^ l+1 W2m+lUJ2m+3 — hm+\Jm+l W 'im+2W2m+i) (61) 

The validity of the above continuity equation (I60j) depends on two assumptions only: 
(i) All Lindblad operators commute with the density H m in the bulk, 2 < m < n — 2 
(second equality sign), and (ii) [H m , H m i] — if \m — m'\ > 2 (third equality sign). 

Using eq. ( 14"T|) of theorem 3 we can now compute NESS expectation values of energy 
density H m and energy current Q m , and also of somewhat simpler spin density 

a z m = -lW2m~lW2m (62) 
and spin current 

S m = CT^(7^ +1 - (7^(7* +1 = -lW2mW2m+2 ~ iW2m-lW 2 m + l (63) 

which are all quadratic in Wj. Note, however, that the spin density satisfies continuity 
equation (d/di)(cr^) = —(S m ) + (S m -i) only in the isotropic case, when = JJ n . 

4-1. Homogeneous transverse Ising chain 

Here we limit our discussion to the spatially homogeneous case J*' y = J x,y , h n = h. We 
shall show that in this case the eigenvalue problem 

Av = (3v (64) 

for (1571) can be most easily and elegantly treated if formulated in terms of an 
abstract inelastic scattering problem in one dimension, with asymptotic solutions 
given in terms of free normal modes for the infinite translationally invariant chain 
v = (. . . ,M^ m_1 ,M^ m ,M^ m+1 , • • -) T , where £ is a complex quasi-momentum (Bloch) 
parameter and u is a 4-dimensional amplitude vector satisfying the condition 

(-Rf r 1 - hK + Ri£ - PU)u = (65) 
and the baths playing the role of inelastic (absorbing) scatterers at the edges of a finite 
lattice. The 'elastic' (Hamiltonian) version of this trick has been used to compute 
temporal correlation functions in kicked Ising chain |3T|. 

The singularity condition for the free mode equation fl65l) results, for a general 
homogeneous XY model, in eight master bands - different values of momenta £ for each 
value of the spectral parameter (rapidity) (3. In order to simplify the discussion - which 
will still get rather involved - we shall in the following restrict ourselves to the transverse 
Ising model J x = J, J y = 0. In this case we find just two master bands with simple 
dispersion relations 

h 2 + J 2 + f3 2 ±u{(3) 



UP) = l hJ — := VW+ .P + P 2 ) 2 - (2hJ) 2 (66) 

but each band is doubly-degenerate, since the corresponding amplitude problem (|65l) 
has two linearly independent solutions 

f-h(h 2 - J 2 + (3 2 ±uj)\ ( \ 





(3(h 2 + J 2 + (3 2 ±u) 




-h(h 2 -J 2 + p 2 ±u) 


\(3(h 2 + J 2 + (3 2 )±lu) ) 



(67) 
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Naively speaking, £_ represents left moving and £ + right moving free modes, each having 
two possible polarizations. Note that = 1. For a general complex (3 we shall choose 
the branch of square root u}(/3) ( 166]) for which |£. 
problem on the left bath in terms of an ansatz 

/ Ml 

^x£ + ipzUz + ipt Mi + ^2^2 
(ipiuj + ip2lh)€- + Of Mi" + i>2 M2 )£+ 



Mi + ^ 2 m 2 )£- + (iPiut + y£)t 



< 1. Let us now write the scattering 
\ 

(68) 



V 



/ 



where u L represents a 4-dimensional vector of left-most eigenvector components, tpi 2 
are the amplitudes of (known) incident free modes, and ip^ 2 are the amplitudes of the 
scattered, outgoing free modes. Plugging the scattering ansatz to the eigenproblem ( 1641) . 
the first two rows of A (1571) result in 6 linearly independent equations for 6 unknowns 
ipi 2 ,u L . Eliminating four variables u L we finally arrive to the non-unitary 2x2 S- matrix 



^2 



(69) 



with 



D n 
°12 

qL 
qL 

°22 



(70) 



r-V 2 (-(r^) 4 + 4(r^) 2 (/3 2 - 3h 2 ) - 16h(hJ 2 + iT L uj)) 
r- 1 /5((r+) 3 + 8iT L _h(3 + AT^h 2 - (3 2 )){2iu) 
r-V((r^) 3 - 8iT L _h(3 + AT^h 2 - f3 2 )){-2\uu) 
r-V 2 (-(r^) 4 + A(T\) 2 {(3 2 - 3h 2 ) - 16h(hJ 2 - iT L uj)) 
T : = (T\ff3 2 + 8p 2 (h 4 + (J 2 + (3 2 ){J 2 + p 2 -cu) + h 2 (2/3 2 - to)) 

- 2(TX)\h A + J 4 + 3/3 4 + J 2 {2(3 2 -uo)- (3 2 oo + h 2 (uj - 2J 2 - A(3 2 )) 

Similarly, one can solve the scattering problem from the right bath with the scattering 
ansatz 

/ \ 

ty+ui + iPtut)t+ + (i>iih + fe" )C 1 ( 71 ) 

1P1 3*1 + V^Mj + V^Mi + V'sTM^ 
V Mr / 

defining the right S-matrix 

^2~ 

Note that since the two directions of free modes f)67p do not have left-right symmetry an 
explicit expression for S R is considerably more complicated than flTOj) and shall not be 
written out here. We shall now show that there exist two qualitatively different types 
of NMM - complex rapidities (3 solving ( 1641) for sufficiently large n. 

First, we shall discuss the so called evanescent normal master modes. These are 
characterized with amplitudes ( 1681) which decay exponentially with the distance from 



V>i + 



(72) 
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Figure 1. Rapidities 0j (black dots) for a transverse Ising chain with J = 1.5, h = 1 
and bath couplings T\ — 1,^ = 0.6, = 1,T^ = 0.3, for three different sizes n= 6 
(upper), n = 30 (middle), and n = 150 (lower panel). Big blue/red dots indicate 
positions of evanescent rapidities (solutions of eq. (|74"l) ) for the left/right bath. 



- say the left - bath, so the other - the right boundary condition becomes physically 
irrelevant in the limit n — > oo. Such solutions if)f 2 = of eq. (|69|) exist exactly when 
the determinant of S-matrix vanishes det S L = 0. Using ([70]) the determinant can be 
written as det S L = ((3/r) 2 p L ((3) wherfjjjj 

p L ((3) = (Tl) s (3 2 - A(T\)\(h 2 -J 2 ) 2 + (2J 2 -Ah 2 )(3 2 + 3/? 4 ) 

- 16(T L + ) 4 (2h 2 (h 2 - J 2 ) 2 - (7h 4 -6h 2 J 2 +2J 4 )[3 2 + A(h 2 -J 2 )f3 4 - 3/3 6 ) 

- Q4(T L + ) 2 (h 4 (h 2 -J 2 ) 2 - 2h 2 J 4 (3 2 - (2h 4 +Ah 2 J 2 -J 4 )f3 4 + 2J 2 f3 6 + (3 8 ) 
+ 25Qh 4 J 4 (3 2 (73) 

% Trivial zero j3 = of course does not represent a physical solution since then the whole S-matrix 
([70f vanishes. 



A general method to solve master equations for quadratic open Fermi systems 



16 



Thus, for sufficiently large spin chains we find at most four NMM whose rapidities are 
given as the roots of 4-th order polynomial in j3 2 

p L (/?evan) = (74) 

that are not simultaneously zeroes of t{(3). Clearly, such NMM asymptotically do not 
depend on the chain size n. In addition, we find evanescent NMM corresponding 
to the right bath simply by replacing by r R in f)74|73l) . In fig. [1] we compare 
evanescent rapidities computed from eq. ( 1741) to numerically calculated spectrum of A, 
at several different sizes n, for a typical case of transverse Ising chain, J = 1.5, h = 1.0, 
strongly coupled to two baths at considerably different temperatures, T\ — 1, T% — 0.6, 
Tf = 1, Tf = 0.3 Note that the same parameter values will be used for numerical 
demonstrations throughout this subsection. 

Second, we shall discuss the other extreme of soft normal master modes with 
rapidities that are closest to the imaginary axis, and thus determining the spectral 
gap of the Liouvillean and relaxation time to NESS. Composing the scattering from the 
two baths with the free propagation along the chain (back and forth) we arrive at the 
general secular equation for the eigenvalue problem f l64l) in terms of a 2 x 2 determinant 

det(£t (n ~ 3) S R S L -I 2 ) =0 (75) 

In the absence of the baths, r ± ' = 0, the solutions of the above problem exist only for 
real quasi-momenta, namely £± = exp(±i , (9),^ e R. For such extended master modes 
the local coupling to the baths can be considered as a small perturbation, thus only 
slightly perturbing the Bloch-like bands /3(e 11? ) = ±ie(i?) with 'energy' 

e(0) = \/h 2 + J 2 - 2\hJ\ cos$ (76) 

The softest NMM, namely the one for which the coupling to the baths is expected to 
be the weakest, should have nearly nodes at the ends of the chain, i.e. d ~ 7r/n, or 
d pa 7i + it /n, and should thus lie near the band edges ±i|/i] ± i| J\ (see fig. [1]). In the 
following we shall focus our calculation on the band edge (3* = i(\h\ + | J|) which, as can 
be checked aposteriori by a straightforward but tedious calculation, always gives smaller 
real part of the rapidity than the lower edge i(\h\ — \J\), and hence really determines 
the gap of the Liouvillean. So we write 

f3 = i(\h\ + \J\) + z (77) 

where z G C is a small parameter, and expand the S-matrices around the band edge 

S L < R = -I 2 + i^Z L ' R V=ii + (78) 



where g := x/^^y, ^ L ' R := (r+' R ) 4 + 4(I^ R ) 2 (4/i 2 + 2\hJ\ + J 2 ) + 16h 2 J 2 and 

z 1 L 1 = 4|/.|(r^) 2 + i6|/ i |(|/ i | + |j|)(|j|-ir L ) 

Z\ 2 = - 2(r L _) 3 - 16T L _\h\{\h\ + \J\) - 8{2h 2 + 2\hJ\ + J 2 ) 
Z\ x = + 2(r L _) 3 - 16r L |/i|(|/i| + \J\) + 8{2h 2 + 2\hJ\ + J 2 ) 
Z\ 2 = 4\h\(Tl) 2 + 16\h\(\h\ + \J\)(\J\ + ir L ) (79) 
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and 

Zf x = (Tl)\2\h\ + \J\) + 4(r R ) 2 (8|/f + 9h 2 \J\ + Mh\J 2 + \J\ 3 ) 

+ 16h 2 \ J\(\J\(3\h\ + 2\J\) - iT K (\h\ + | J|)) 
Z£ = - 2(r*) 3 - \QY K \h\(\h\ + | J|) - 8(2/i 2 + 2|/iJ| + J 2 ) 
Z* = + 2(r R ) 3 - 16r R |/i|(|/i| + | J|) + 8(2h 2 + 2|/iJ| + J 2 ) 
Zf 2 = (Tl)\2\h\ + \J\)+ 4(r R ) 2 (8|/f + 9h 2 \J\ + A\h\J 2 + | J| 3 ) 

+ 16/i 2 | J|(| J|(3|/i| + 2\J\) + ir R (|/i| + \J\)) (80) 

Next we expand £ + ( 1661) in z, yielding 

i + = -l-g- l ^z + 0(\z\) (81) 

and so the free propagator in ( 1751) can be written as 

e + in ~ 3) = exp(2^- 1 v / ^ii) + 0(\z\). (82) 

In eqs. f J7H|81ll82l) the branch cut along the negative real axis has been chosen for \/—iz. 
Since the product of S-matrices in (1751) is near identity, the free propagator should be 
near one as well, hence 2ng~ 1 \ / — iz should be near 27ri. Let us define Zq by setting 
2ng^ 1 y / — iz = 27ri, so 

Zo = -m 2 g 2 n - 2 (83) 

and write z = z (l + y) where \y\ <C 1 is another small complex parameter. However, 
since z is purely imaginary, we need to compute a small but non-vanishing y which 
will, in the leading order in n, solve ( 1751) since the real part of the soft mode's rapidity 
is determined as 

Re (3 = Re z y = n 2 g 2 n~ 2 Im y (84) 



Now, writing y/—iz = ^—iz Q ^l + y = mgn^ 1 ^ + y/2 — y 2 /8) + 0(y 3 ) in 
(I78|82l) . plugging all that to eq. (1751) and computing to order (D(n~ 2 ), noting that 
= 0(n~ 2 ), we arrive to a simple quadratic equation for y, whose solution, plugged 
to (184"]) . gives the final result, namely the sectral gap of Liouvillean A = 2Re/3 
_ (2M) 2 Ai 3 4 
A " (\h\ + \J\)* A™ +U[H } (85j 
Ai := 64(1^ + Tl)h 2 J 2 (2h 2 + 2\hJ\ + J 2 ) 
+ i6((r^) 3 + {Tlf)h 2 J 2 

+ i6rjr;(ri + r R )(2/i 2 + 2|/ij| + J 2 )(4h 2 + 2\hj\ + j 2 ) 

+ 4r L r R ((r L _) 3 + (r^fr* + r^r*) 2 + (r*) 3 )(2/> 2 + 2\hj\ + j 2 ) 

+ (r^r;) 3 (r^ + r;) 

A 2 := ((r^) 4 + A(T h + ) 2 (Ah 2 + 2\hJ\ + J 2 ) + lQh 2 J 2 ) 
x ((r R ) 4 + 4(r*) 2 (4/i 2 + 2\hJ\ + J 2 ) + 16h 2 J 2 ) 

In fig. (j2J) we compare this analytical result to exact numerical calculations of the 
eigenvalue of A with minimal real part, confirming both, its precise numerical value 
and that the relative scaling of the next order correction is indeed (9(ri _1 ). 
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Figure 2. Spectral gap A times a third power of the chain length n for a transverse 
Ising chain with J = 1.5, h = 1 and bath couplings 1^ = 1, 1% = 0.6, Tf = 1, Tf = 0.3. 
Thin horizontal line indicates the theoretical asymptotic value (|85[) . In the inset we 
show deviation from asymptotic constant value of An 3 in log-log scale and demonstrate 
that it decays as oc n^ 1 (thin line). 
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Figure 3. Complete spectrum of 2 12 complex eigenvalues of Liouvillean for a 
transverse Ising chain with n = 6 spins and J — 1.5, h — 1 and bath couplings 
T\ = 1,1% = 0.6, rf = l,Tf = 0.3 (the case of the upper panel of fig. [p). 



Note that, interestingly, both main analytical results of this subsection, namely 
evanescent and soft mode rapidities do not depend on r^ ,R . Physically speaking, they 
only depend on the effective strengths of the bath couplings and not on the temperatures. 

We end this subsection by presenting some further numerical results on heat 
transport in the open transverse Ising chain in the Lindblad form. In fig. [3] we 
demonstrate expression (1421) for constructing the full spectrum of the Liouvillean in 
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Figure 4. Energy current (upper/blue points), and average spin current (lower/red 
points), versus chain length n for a transverse Ising chain with J = 1.5, h = 1 and bath 
couplings = 1,1% = 0.6, Tf = 1, Tf = 0.3. 



0.0 




Figure 5. Energy density profile (lower, blue points), and spin density profile 
(upper, red points), for a transverse Ising chain of n — 80 spins with J = 1.5, h = 1 
and bath couplings T\ = l,r£ = 0.6, Tf = l,Tf = 0.3. The insets display 
logarithm of the difference to the bulk values 8H m := \(H m ) — -ffbuikl (blue points), 
<5<7^ := — crg ulk | (red points) in comparison with ±(41og£_)rn + const with 

quasi-momentum £_ = 0.584692 corresponding (]66[) to the leading evanescent rapidity 
/3 ovan = 0.438739 (full lines). 



terms of a set of rapidities, for a short chain. In fig. H] we demonstrate eq. (HTj) of 
Theorem 3 by computing the energy current Q m (ED), and the average spin current 
S = Y^m=x (|63|) in NESS of a typical transverse Ising chain. Numerical results 
give a clear indication of ballistic transport (Q) = O(n ), (S) = O(n ), however its 
rigorous proof and analytical calculation of the currents would require full control over 
the complete set of NMM which is at present not available. In fig. Owe plot the energy 
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Figure 6. Average Liouvillean spectral gap (A) versus the chain length n for 
disordered XY models: (i) Jf n = 0.5, = 0, h m € [1,2] (transverse Ising with 
field disorder, blue points), (ii) € [0.5,2], J-jJ, = 0, h m = 1 (transverse Ising 
with interaction disorder, red points), (iii) e [0.5,1], € [0.5, l],/i m = 1 (XY 
with interaction disorder, golden points), all for bath couplings T\ = 1,^ = 0.6, 
r^ 1 = 1, rJ? = 0.3. Full lines indicate exponential fits to right halves of data. Averaging 
is performed over 2000 disorder realizations. 

density (1591) and spin density (1621) profiles in NESS. Again, we note flat profiles in the 
bulk of the chain, m,n — m ^> 1, with exponential falloff due to adjustment to the 
non-equilibrium bath values. Since the densities can be written, by means of (j4"T|) . 
as 4— point functions in NMM components, the leading falloff exponents of the profile 
\(H m ) — i?bulk| ~ |£-| 4m is given by the quasi-momentum ( 1661) corresponding to the 
maximal evanescent rapidity /3 evan (1741) . 

4-2. Disordered XY chain 

In this subsection we treat the opposite extreme, a disordered XY chain (|52|) where three 
sets of physical parameters are chosen as random uncorrelated variables from uniform 
distributions on the intervals, J* e [J^, J^J, J y m G [JL^maJ. G [h min, "'max I • 
Clearly, the eigenvalue problem (|64p for the matrix (|57|) then becomes equivalent to the 
Anderson tight-binding problem in one dimension for a quantum particle with a 4— level 
internal degree of freedom. We do not pursue any theoretical analysis of this problem 
here, but merely state that numerical investigations indicate existence of exponential 
localization of all eigenvectors (or normal master modes) for disorder of any strength in 
anyone of system's parameters. 

With the picture of localization of NMM in mind, the effect of the couplings to the 
heat baths at the chain's ends on quantum transport can be predicted by theoretical 
arguments (see [32] for a review of related studies): The spectral gap of the Liouvillean 
should be exponentially small A ~ exp(— n/£) where I is the localization length of NMM 




A general method to solve master equations for quadratic open Fermi systems 



21 




10 20 30 40 50 



n 

Figure 7. The scaling of the energy current (Q m ) with chain length n of the disordered 
XY model in the same regimes/parameters/plot styles as in fig. [5] 
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Figure 8. Scaled energy density profile of interaction disordered XY chains (case 
(iii) of fig. [5]) for three chain sizes: n = 20 (blue points), n = 40 (red points), n = 60 
(golden points). Averages over 50000 disorder realizations have been performed. 

which is expected to be proportional to the square of inverse disorder strength. This 
is demonstrated in figj6j If all NMM are exponentially localized, the currents should 
decrease with the chain size n faster than any power, perhaps exponentially, and the 
system should behave as an ideal insulator (at all temperatures). This is demonstrated 
by straightforward numerical calculations of the heat current (1611) in fig. [71 In the final 
figure M we plot the energy density profile (H m ) (159|) in a typical case of disordered XY 
chain, versus a scaled spatial coordinate (m — l)/(n — 1) G [0, 1], for several different 
chain lengths n, and demonstrate sharping up of energy density profiles with increasing 
n, which is again indicating insulating behaviour. 
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5. Discussion and conclusions 

The main result of the paper is a general method of explicit solution of master equations 
describing dynamics of open quantum system, under the condition that the system's 
Hamiltonian is quadratic and all Lindblad operators are linear in canonical fermionic 
operators (which can either represent real physical fermions or any abstract 2-level 
quantum systems (qubits) thru the Jordan- Wigner transformation). Using a novel 
concept of Fock space of physical operators (or density operators of physical states), and 
the adjoint structure of canonical creation and annihilation maps over this space, the 
problem can be treated in terms of a non-Hamiltonian generalization of the method of 
Lieb, Schultz and Mattis [18] lifted to an operator space. We have explicitly constructed 
a non-canonical analog of Bogoliubov transformation of the quantum Liouville map to 
normal master modes. Related ideas in the Hamiltonian context have been used by the 
author [31], [331 El] in order to approach the problem of real time dynamics and ergodic 
properties of isolated interacting many-body quantum systems. 

As an illustration of the method we have solved far from equilibrium quantum 
heat and spin transport problem in Heisenberg XY spin 1/2 chains which are coupled 
to canonical heat baths only at the two ends. Irrespectively of the strength of the 
coupling to the baths and their temperatures, we have shown a ballistic transport in 
the spatially homogeneous (non-disordered) case, and an ideally insulating behaviour in 
the disordered case associated to localization of normal master modes of the quantum 
Liouville operator. In this context the method can be considered as a simple alternative 
to the solution of quantum Langevin equations |24j . 

However, the method should easily be applicable to variety of other physical 
situations, for example if all fermions are coupled to the baths one could make a 
solvable model of genuine quantum diffusion, a many-body generalization of the tight- 
binding model [35J. We also expect the method to be applicable to the Redfield type 
of master equations (see e.g. [35]) - which do not conserve positivity for a short initial 
(slippage) time interval - provided only the system part of the Hamiltonian is quadratic 
and system-bath couplings are linear in fermionic variables. Furthermore, extension of 
the method to open many-boson systems should be straightforward, simply by replacing 
anticommutators by commutators throughout the exposition of section [21 

Treating density operators of NESS as elements of a Hilbert space of operators one 
may also extend the concept of entanglement entropy, with respect to a bipartition of a 
system of many fermions [36], to NESS which can in our approach be viewed as a kind of 
ground state of the Liouvillean. Saturation of such operator space entanglement entropy 
[31] (which is suggested by numerical experiments [37]) indicates efficient simulability 
of NESS by elaborate numerical methods such as density matrix renormalization group 
[38] . perhaps even for more general, non-solvable quantum systems. 

As last we mention a more ambitious extension of the present work: Namely we 
propose to explore a question, whether more involved algebraic methods of solution of 
interacting many-body quantum systems, like e.g. Bethe Ansatz or quantum inverse 
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scattering [IB] , could be generalized to open quantum systems, e.g. by means of the 
proposed concept of Fock space of operators. Could one discuss completely integrable 
open quantum systems which go beyond quadratic Liouvilleans? 
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